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COMPUTATION OF ORBITS USING TOTAL ENERGY 


1.0 SUMMARY 

The computation of orbits can be done more efficiently by the use of any of several 
new formulations (Reference 1, 2, 3, 4, 5) of the perturbed two body problem which 
consider the total energy of the orbital system as one of the dependent variables. The 
total energy is the osculating two body energy plus the potential energy due to 
perturbing masses. The use of the total energy as the dependent variable instead of 
the two body energy is a relatively new idea (Reference 1). The advantage of using 
total energy arises from the fact that the more perturbing potential energy that is 
accounted for in the total energy variable, the more nearly constant is the total 
energy. In fact, except for dissipative forces such as drag, the only reason for the 
total energy not being constant is the rotation or revolution of the perturbing mass. 
This near constancy of the total energy has the effect of inhibiting error growth 
during numerical solution (Reference 1). This paper will present the results of an 
application of total energy formulation (Reference 2) to the problem of the precise 
computation of orbits. 

2.0 INTRODUCTION 

The differential equation of motion of the perturbed two-body problem can be 
expressed as, 


r 



r = F = P - 


dV 

dr 


(D 


where r is the position vector of one of the bodies relative to the other. The perturba- 
tions, those derivable from a potential dV/d r, as well as other forces P, are included in 
the total perturbation F. 


The total energy element formulations (References 1, 2, 3, 4, 5) of the perturbed two- 
body problem are developed such that dV/dt is used as well as dV/dr. The perturba- 
tions are split into those derivable from a potential and those which are included in 
the perturbation P. The perturbation P normally includes non-conservative pertur- 
bations, but it can also include perturbations derivable from a potential when 
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convenient. For the total energy formulations, the right-hand sides of all differential 
equations, except that for the total energy, include the perturbation factors P and 
dV/dr. The total energy differential equation has the form 


h 



+ 


av 

at 


where h is the total energy, 

1 • • u 

h = - r.r - — + V(r,t) 

2 r 


( 2 ) 


(3) 


Equation (2) is derived by taking the time derivative of equation (3) and sub- 

• • 

stituting equation (1) into this result to eliminate r. Note that this differential 
equation includes the perturbations P and dV/dt, but does not include dV/dr. 


There are three options available in the total energy formulations depending upon 
the way in which the perturbations derivable from a potential are used in the 
differential equations. These options are categorized as follows: 

(A) The entire perturbing potential is considered with its effect including dV/dr and 
dV/dt. This is the option which is developed and discussed in this paper. 

(B) The perturbing potential can be portioned, including some of the perturbation 
in dV/dr and some in P. This has been the approach most often used when the 
geopotential is the perturbation. The zonal terms have been included in dV/dr, 
while the explicitly time dependent terms (the tesseral and sectorial terms) 
have been included in P. This approach has been used in order to avoid the 
computation of dV/dt. The potential used is that of the zonal terms only. 

(C) The perturbing potential is not considered at all. The dV/dr are included in the 
perturbation P. The potential is set to zero. 
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It must be emphasized that all three options are correct. The advantage that any of 
the options has over the others is numerical accuracy and speed in computation. The 
advantage of option (B) over option (C) in accuracy and speed is considerable and is 
discussed at length in References 1, 2, and 4.* 

In order to properly implement the differential equation (2) for the total energy as 
discussed in option (A), the partial derivative dV/dt must be computed. This report 
will derive a simple formula for this computation. This formula will be developed for 
the general case of any perturbation derivable from a potential. Then the particular 
case of a geopotential perturbation acting on an Earth satellite will be used as an 
example to show the advantage of using dV/dt in the computation. 

3.0 DEVELOPMENT OF aV/at 

Let r be the position vector in an inertial system and let r G be the same position 
vector in a system rotating with angular velocity co with respect to the inertial 
system. Then, 

r = Ig (4) 

The velocity vectors are related by, 

• • 

r = r G + co x r G (5) 

In the inertial system, the potential function is expressed as an explicit function of 
time, 

V = V(r,t) (6) 


* For these formulations a slightly different energy parameter a 0 , where 
h = — 2a 0 , is used and a new independent variable called fictitious time is 
introduced. With these changes, equation (2) becomes 



where ( )' = d( )/ds and s is the independent variable such that dt/ds = r. 
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having the total derivative 


dv = av . ; + jv 

dt ar at 


(7) 


In a properly chosen rotating system, the same potential function can be expressed 
as a function of position only, 

V = V(r G ) (8) 


having the total derivative, 

iV = JV . * G (9) 

dt ar 

since in the rotating system the potential has no explicit dependence on time, 

av(r G )/at — o. 


Using equations (4) and (5), equation (7) becomes 


d V 
dt 


— • ( r G + co x r G ) + 

ar G At 


( 10 ) 


Comparing equations (9) and (10), we obtain 

• CO x r G (11) 

at ar G 

Note that to this point, we have not considered any particular potential function. 
The result, equation (11), can be applied under proper conditions to the case where V 
represents the perturbing geopotential function or to the case of a lunar, solar, or 
planetary perturbation on a satellite. 


4.0 APPLICATION TO THE GEOPOTENTIAL 

We now consider the case of the perturbing geopotential which can be divided into 
two parts, 

V(r,t) = V z (r) + V T (r,t) (12) 


where r is expressed in an inertial system having one axis normal to the Earth 
equatorial plane and the other two orthogonal axes in the equatorial plane. The 
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portion of the perturbing geopotential Vz(r) arising from the zonal terms has no 
explicit time dependence. The portion of the perturbing geopotential V T (r,t) arising 
from the sectorial and tesseral terms are explicitly dependent upon time. 

For the case of a perturbing geopotential equation (11) can be reduced further. 
Define the rotating system X G , Y g , Z g such that Z G is in a direction normal to the 
Earth equatorial plane and the X G and Y G axes lie in the Earth equatorial plane and 
are fixed in the Earth. The zonal portion of the perturbing geopotential is 

v z=- -S- 2 C n ,o^-^-) P n (Zc/r) 
n — 2 


where, C n ,o are the zonal coefficients 

a e is the equatorial radius of the Earth 

P n is the n th degree Legendre polynomial which is a function of Zc/r. 


Now, r G — i G X G + j G Y G + k(jZ G , 



where i G , j G , k G are unit vectors along the X G , Y G , Z G axes. The partial 
derivative dVz/dr G has the form, 


av z 

dr G 


= f^r.Zolro + f 2 (r,Z G )k G 


Note that since co = cok Gj 

. co x r G = 0. 

ar G 


15b 



Thus , the zonal part V z of the perturbing geopotential does not contribute to dV/dt. 
Equation (11) becomes 


= - • wxr G , 

dt dr G 

Consider the first two options given in Section 2.0 for the formulation of the 
differential equation (2) for the total energy. 

4.1 Option A - All zonal, sectorial, and tesseral terms of the perturbing 
geopotential are included in the potential and hence in the total energy. 

Let V = V z (r) + V T (r) 


then 


also P = 0 

and from equation (13), we compute dV/dt. 

Further, since r = r G , we can express equation (2) as 


h A = 


. co x r 


4.2 Option B - Only the zonal terms of the perturbing geopotential are included in 
the potential and hence in the total energy. 

Let V = V z (r) 


then 
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The sectorial and tesseral terms are considered to be in the perturbation P, 


p _ _ dV T 

dr 

and the total energy differential equation (2) becomes, 

* * dVr " 

h B = r • P = - — - . r (15) 

dr 


4.3 Comparison of Options A and B : 

Both equations (14) and (15) depend directly upon the factor dVr/dr, which is a small 

term depending only upon the sectorial and tesseral terms. But we also observe that 
• • • 
for Option B, h B is proportional to the inertial velocity, r, whereas for Option A, h A 

is proportional to the component (co x r) of the inertial velocity which arises from 

the rotation of the axes fixed in the Earth. 

For near Earth satellite orbits, 




• 

co x r 

< 

r 





and also 


(16) 


• 


• 

h A 

< 

hB 


In fact, if the Earth were not rotating (co = 0), then h A would be zero. For satellite 
orbits which are at large distances from the Earth, the inequality (16) does not 
always hold. However, at large distances from the Earth, the perturbing 
geopotential is not as significant as perturbations due to the Sun or Moon. The 
global region for which the inequality (16) holds is complicated and depends upon the 
semi-major axis, the eccentricity, and the true anomaly (or angular position) of the 
satellite trajectory. 
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For near circular orbits, it can be shown that the ratio 


\ 1/2 

1 - (Z/r) 2 j 

where n is the mean motion of the satellite and the factor (l - (Z / r ) 2 ) is always 
less than unity. For orbits within the geosynchronous distance, the inequality 
(16) holds since 

n > co 

For near Earth orbits, 

co _ 1 

n 16 

and so 

h A /h B < — 

16 

The formulations of the perturbed two-body problem discussed in References 1, 2, 
and 4 are in effect perturbed harmonic oscillators having frequencies which are 
dependent upon the total energy. The use of the full geopotential as shown in 

Option A in the computation of the total energy causes h A to be small. Thus, h A is 
nearly constant and the resulting frequency of the perturbed oscillator equations is 
nearly constant. Options A and B as well as Option C are also compared in Table I. 

5.0 NUMERICAL RESULTS 


h A / h B = 


The numerical effect of using the full geopotential as in Option A is shown in Figure 
1. A near circular orbit was propagated for ten days first using Option A and then 
using Option B. This computation was done using the KSUR12 total energy formu- 
lation (Reference 2) and the RK4(5) variable step numerical integrator (Reference 6). 
The geopotential model used was the complete GEM-L2 (Reference 7). The results of 
these computations were compared to a reference trajectory computed with very high 
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precision as given in Reference 8 and originally provided in Reference 9. Figure 1 
shows the RSS of the position vector of Options A and B, with each compared to the 
reference. 

Option B (using only the zonal part of the geopotential in the total energy) required 
an average of 59.4 variable steps per revolution with a maximum error of 25 meters. 
Option A (using the full geopotential in the total energy) required an average of 45.2 
variable steps per revolution with a maximum error of about 8 meters. The two 
options are also compared on Figure 1 using 30 fixed steps per revolution. Option B 
showed a rapidly growing error reaching 25 meters after 4 days and still diverging. 
Option A reached a maximum error of about 15 meters after TO days. 
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Table I. Comparison of Options A, B f and C for 
perturbing geopotential. 



Option A 

Option B 

Option C 

Perturbing 

Geopotential (V) 

Vz + Vt 

V z 

0 

Total Energy (h) 
equation (3) 

1 * • U 

— r.r + V z + V T 

2 r 

1 ; ; ii , v 

— r 1- Vz 

2 r 

1 * * pi 

— r.r - 

2 r 

u • ( dY \ 

Perturbation l 1 

v ar / 

dVz 3Vx 

dr dr 

av z 

ar 

0 

Perturbation (P) 

0 

dV T 

dr 

1 av z av T | 

-[* + ar J 

Derivative of # 
total energy (h) 
equation (2) 

3Vt 

- • co x r 

dr 

av T 

- • r 

dr 

( 3V Z 3V t \ • 

\ dr + ar J 
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FIXED = fixed step RK5 

VAR = variable step with RK4(5) 

SPR = steps per orbital revolution 
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Figure 1 . Comparison of RSS errors for near Earth trajectories computed using only the 
zonal terms of the perturbing geopotential in the total energy (Option B) and 
using the full perturbing geopotential in the total energy (Option A). 
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